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Abstract 

We construct a model of speciation from evolution in an ecosystem consisting of a 
limited amount of energy recources. The species posses genetic information, which is 
inherited according to the rules of the Penna model of genetic evolution. The increase 
in number of the individuals of each species depends on the quality of their genotypes 
and the available energy resources. The decrease in number of the individuals results 
from the genetic death or reaching the maximum age by the individual. The amount of 
energy resources is represented by a solution of the differential logistic equation, where 
the growth rate of the amount of the energy resources has been modified to include 
the number of individuals from all species in the ecosystem under consideration. 

The fluctuating surrounding is modelled with the help of the function V(x,t) — 
\x* + ^b(t)x 2 , where x is representing phenotype and the coefficient b(t) shows the 
cos(ui) time dependence. The closer the value x of an individual to the minimum of 
V(x, t) the better adapted its genotype to the surrounding. We observed that the life 
span of the organisms strongly depends on the value of the frequency ui. It becomes the 
shorter the often are the changes of the surrounding. However, there is a tendency that 
the species which have a higher value an of the reproduction age win the competition 
with the other species. 

Another observation is that small evolutionary changes of the inherited genetic 
information lead to spontaneous bursts of the evolutionary activity when many new 
species may appear in a short period. 



1 Introduction 

There have appeared many papers and books on the dynamics in ecological 
systems. Many of them start from the Lotka-Volterra model ^] which de- 
scribes populations in competition. Usually the populations are represented by 
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the various types of the predator-prey systems. They may exhibit many inter- 
esting features such as chaos (e.g. the recent papers on the topic M - and 
phase transitions (e.g. [ p[ |7, The possible existence of chaos became evident 
since the work of May [pTlffl. However, usually the studies concerning chaos 
in biological populations do not contain the discussion of the role of the inher- 
ited genetic information in it. The discussion of a predator-prey model with 
genetics has been started by Ray et al. pL who showed that the system passes 
from the oscillatory solution of the Lotka-Volterra equations into a steady-state 
regime, which exhibits some features of self-organized criticality (SOC). Our 
study [|| on the topic was the Lotka-Volterra dynamics of two competing pop- 
ulations, prey and predator, with the genetic information inherited according 
to the Penna model [|l2) of genetic evolution and we showed that during time 
evolution, the populations can experience a series of dynamical phase transi- 
tions which are connected with the different types of the dominant phenotypes 
present in the populations. Evolution is understoood as the interplay of the 
two processes: mutation, in which the DNA of the organisms experiences small 
chemical changes, and selection, through which the better adapted organisms 
have more offsprings than the others. The problem of speciation from evolution 
has been studied recently by McManus et al. jl3| in terms of a microscopic 
model. They confirmed that the mutation and selection are sufficient for the 
appearance of the speciation. We followed their result and in this study we 
consider a closed ecosystem with a variable number of species competing for 
the same energy resources. The total energy of the ecosystem cannot exceed 
the value f2 and every individual costs a respective amount of energy units. We 
show that the small evolutionary changes of the inherited genetic information 
result in the bursts of evolutionary activity during which new species appear. 
Some of the species become better adapted to the fluctuating surrounding and 
they win the competition for the energy resources. 

Our model belongs to the class of Lotka-Volterra systems describing one 
prey (represented by a self-regenerating energy resources of the ecosystem) and 
a variable number of the predators competing for the same prey. 



2 Evolution of energy resources 

All species in the ecosystem under consideration use the same energy resources. 
In the model, the number -/Vg(i) of the energy units available for the species 
satisfies the differential equation 

dN E (t) (A( N(t) Ne® m 

with the initial condition 

^N E (t )=a, (2) 

where N(t) represents the total number of all individuals in the ecosystem (from 
all species), Se is the regeneration rate coefficient of the energy resources, 7^ 
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represents relative decrease of e_g caused by the living organisms. The equation 
means that regeneration of the energy resources in the ecosystem takes some 
time and its speed depends on the number of the living organisms. The above 
equation also ensures that the size of the species cannot be too big and the 
number of the species coexisting in the same ecosystem is limited. Otherwise 
they would exhaust all the energy resources of the ecosystem necessary for life 
processes. 

In the case when N(t) = 0, i.e., if there are no living organisms in the 
ecosystem, the Eq.[j] reduces to the well known logistic differential equation 
with the following analytical solution Jl3 

1 a 

n NE ^ = a+(l-a)exp - £E {t-t y (3) 

and the initial condition Eq.|[ 

We adapt the above solution into our model (Eq.|l|). To this end we assume 
that the individuals from all species in the ecosystem under consideration may 
reproduce only at discrete time t — 0,1,2,3,... (otherwise N(t) = const), 
whereas NE(t) remains a continuous function of t between these discrete time 
values. Say, if there is N(to) individuals at the discrete time value t = to then 
the value N(t) remains constant (N(i) — N(to)) in the whole time interval 
[to, to + 1). The analytical solution of Eqjl] in this time interval is the following 



h NE ® = ~~7^ ^ 7^ 7\ for 1 < h ' (4) 

" a+(l-a)exp-E E {l-'yE-5 J -){t-to) 

where the initial condition is represented by Eqj^ and t\ = to + 1. At time 
t = t\ the species reproduce themselves and a new value, N(t±), becomes the 
initial condition for the next time interval, [ti,ti + 1). The numbers N(t) 
(t = 0, 1, 2,3,.. .) result from the computer simulation of mutation and selection 
applied to the species in the ecosystem according to the Penna model jl2] of 
genetic evolution. The Penna model of evolution represents evolution of bit- 
strings (genotypes), where the different bit-strings replicate with some rates 
and they mutate. Hence, in our model we have two types of dynamics, the 
continuous one for the number A^(f) of energy units and the discrete one for 
genetic evolution of bit-strings. 



3 Species evolution 

We restrict ourselves to diploid organisms and we follow the biological species 
concept that the individuals belonging to different species cannot reproduce 
themselves, i.e., they represent genetically isolated groups. The populations 
of each species are characterized by genotype, phenotype and sex. Once the 
individuals are diploid organisms there are two copies of each gene (alleles) 
in their genome - one member of each pair is contributed by each parent. In 
our case, the genotype is determined by 2L alleles (we have chosen L = 16 in 
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computer simulations) located in two chromosomes. We agreed that the first 



L' sites in the chromosomes represent the housekeeping genes 15 , i.e. the 



genes which are necessary during the whole life of every organism. We assumed 
that there also exist L — L' additional "death genes", which are switched on 
at a specific age a = 1,2, .... L — L' of living individual. The idea of the 
chronological genes has been borrowed from the Penna model |l^, [l7[] of 
biological ageing. The term "death gene" has been introduced by Cebrat |L8[ 
who discussed the biological meaning of the genes which are chronologically 
switched on in the Penna model. The maximum age of individuals is set to 
a = L — L'. It is the same for all species. All species have also the same number 
L' of the housekeeping genes and the same number L — L' of the chronological 
genes. However, they differ in the reproduction age or, i.e., the age at which 
the individual can produce the offsprings. The individuals can die earlier due 
to inherited defective genes. We assume, that it is always the case, when an 
individual reaches the age a and in its history until the age a there have appeared 
three inactive genes in the genotype (inactive gene means two inactive alleles). 

We make a simplified assumption, that all organisms fulfill the same life 
functions and each function of an organism (does not matter what species) has 
been coded with a bit-string consisting of 16 bits generated with the help of a 
computer random number generator. We have decided on L different functions 
and the corresponding bit-strings represent the patterns for the genes. The 
genes of all species are represented by bit-strings which differ from the respec- 
tive pattern by a Hamming distance H < 2. Otherwise the bit-strings do not 
represent the genes. In order to distinguish the species, we have introduced the 
concept of an ideal predecessor, called Eve, who uniquely determines all indi- 
viduals belonging to the particular species. Namely, the individuals have genes 
which may differ from the genes of Eve only by a Hamming distance H < 1. 
Genes, which differ from the genes of Eve by a Hamming distance H > 1 and, 
simultaneously, which differ from the pattern by a Hamming distance H < 2 
represent mutated genes. They are potential candidates to contribute to a new 
species but they are considered as the inactive genes for the species under con- 
sideration. If there happens another individual with the mutant gene in the 
same locus then these two mutants can mate (different sex is necessary as well 
as the age a > an) and they can produce offerings. In the latter case the new 
species is created with Eve who represents Eve of the old species except for the 
mutant genes. In our computer simulations the new species ususally are extinct 
due to the mechanism of genetic drift. However, after long periods of 'stasis' 
there happen bursts of the new species which are able to live for a few thousands 
of generations or even more. They also can adapt better to the surrounding and 
they can dominate other species. The small changes of the inherited genetic 
information are realized through the point mutations. 

A point mutation changes a single bit within the 16-bit-string representation 
of a gen and according to the above assumptions the gene affected by a mutation 
may pass to one of three states: 5 = 1 (gene specific for the species), 5 = 2 
(mutant gene, potential candidate of a new species), 5 = (defected gene). In 
the model, the bit-string representing a defected gene (5 = 0) can be mutated 
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Figure 1: Phenotype-surrounding interaction function V(x,i) 
for two values of b(t) (0.5 and -0.4). 
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with probability p = 0.1. Hence, there is still a possibility for back mutations. 
Genes are mutated only at the stage of the zygote creation - one mutation per 
individual. 

The life cycle of diploids needs an intermediate stage when one of the two 
copies of each gene is passed from the parent to a haploid gamete. Next, the 
two gametes produced by parents of different sex unite to form a zygote. 

In the model, phenotype is defined as a fractional representation of the 16- 
bit-string genes, where each gene is translated uniquely into a fractional number 
x from the interval [—1, 1]. The phenotype is coupled with the surrounding with 
the help of the function V(x, t) as follows: 



V(x,t) = ^+^b(t)x 2 (5) 

where the fluctuations of the surrounding are represented by the coefficient b(t). 
The function V(x, t) has one minimum, x — 0, for b > and two minima, 
x = ±y/\b\, for b < as in Figj|. We use them to calculate gene quality, q, in 
the variable surrounding 



q = e -<Y(*.t)-Vn, 



96 [0,1], 



(6) 



where the parameter T has been introduced to control selection. The smaller 
value of T the stronger selection. Next, we calculate the fitness Q of the indi- 
viduals, say at age a, to the surrounding with the help of the average 



L'+a 



Q 
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L' + a 



(7) 



with 
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q t = max {qVsisV, l),g? 2) 5(S| 2) , 1)}, 



(8) 



where the indices (1) and (2) denote, respectively, the first allele and the second 
allele at locus i = 1, 2, . . . , V + a and Si = 0, 1, 2 represent their states. The 
symbols SQ denote Kronecker delta. The inactive alleles and mutant genes do 
not contribute to Q and always Q £ [0, 1]. 

We decided to determine the amount of possible offsprings by projecting 
the fitness of the parents, Qf (female) and Qm (male), onto the number of 
produced zygotes 



where at least one zygote is produced and the maximum number of zygotes 
is equal to 10. In the computer simulations, the zygotes are mutated, after 
they are created, and they can be eliminated if at least one housekeeping gene 
becomes inactive. The individuals may reproduce themselves after they reach 
the age a — cir. The genotypes of the parents who are better adapted to 
the surrounding generate more offsprings. The mutants, if they happen (they 
posses genes with the Hamming distance H > 1 from Eve) cannot mate with 
the non- mutants. 

We assume that the sex of a diploid individual is determined with the help of 
a random number generator at the moment when it is born and it is unchanged 
during its life. 

4 Computer algorithm 

We investigate the species with respect to speciation. In most computer simula- 
tions we observed evolving single species and analyzed the offsprings of the new 
mutant species originating from the old one. The secondary order speciation, 
i.e., the speciation taken place from the mutant species, has not been consid- 
ered. However, we investigated the case when initially there are a few species 
in the ecosystem and the results qualitatively were the same as for single initial 
species. 

In the simulation, it is very important to prepare the genetic information in 
a proper way. It is obvious that in the evolution process some genes are not nec- 
essary during the whole life, e.g., they may become important only near the end 
of life, and it is possible that they could be defective since the individual under 
consideration was born. Therefore, first we prepare the initial species in the 
time independent surrounding for a few thousands generations until the inher- 
ited genetic information is represented by a steady state flow between suceeding 
generations. Then the distribution of inactive genes becomes time independent. 
Only after that we switch on the fluctuations of the surrounding, which in our 
case arc represented by the coefficient b(t) = — \ — cos{tot), and all the initially 
prepared species are put together in the ecosystem. Hence, the computer algo- 
rithm consists of two blocks: preparation and evolution in variable surrounding. 



N Z y g ote = max{l, 10 x min(<2_F, Q M )} 



(9) 
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INITIAL SPECIES PREPARATION: 



(1) generate L bit-strings (16 bits) representing life functions of the or- 
ganism and which are the patterns for genes 

(2) generate N s predecessors (Eve) of the initial species (Hamming dis- 
tance of each gene from the respective pattern, H <= 2, and then 
S 1 = 1,S 2 = 1,...S L = 1) 

(3) for each Eve, representing the species 1, 2, . . . , N s construct Nm males 
and Np females (Nm = Np) who have each gene at distance H < 1 
from the corresponding gene of Eve. 

(4) evolve separately each species 1, 2, . . . , N 8 as follows: 

(i) increase the age of all individuals by 1 

(ii) remove from the population all individuals who should die be- 
cause they have exceeded the maximum age (L — L') or during 
their life until now they have collected three inactive "chrono- 
logical" genes. 

(iii) [optionally] the Verhulst factor is applied, i.e. every individ- 
ual survives with the same probability (1 — (Nm + Np)/N max ), 
where N max is the maximum allowed number of individuals in 
the population under consideration. 

(iV) select at random a female and a male for whom a > and let 
them produce N zygotes according to Eq|j. 

(V) each zygote is mutated. We assume that the housekeeping genes 
cannot be inactive. Therefore, at least one defected housekeeping 
gene kills zygote. The total numer of the zygotes in the popu- 
lation which survive cannot exceed the number s PO p(Ne/^)Nr, 
where Nr is the total number of pairs of individuals, female 
and male, whose age a > clr. The parameter e pop controls the 
allowed number of zygotes. 

(Vi) Decrease the total number of energy units by the number N(t) 
of all individuals (parents and their new-born offsprings), i.e., 
N E (t) = N E (t-l)-N(t) 

(Vii) regenerate the ecosystem energy according to Eq.|| 
(Viii) goto (i) unless the stop criterion is fulfilled 

SPECIES IN VARIABLE SURROUNDING: 

(5) all N s species are put together in the ecosystem. 

(6) t = t + l 

(7) calculate b(t) representing surrounding fluctuation 

(8) check all individuals from the old species with respect to mutants 
(who have, at least at one locus, k — 1, . . . , L, two bit-strings rep- 
resenting the pair of the alleles at state SW = 2 and =2). If 
there is in the population an individual representing a mutant then 
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create Eve, the pattern of new species, and look for other mutants in 
the same species consistent with the created Eve. Determine a new 
value of the reproduction age an E [1,L — L'\ for the species with 
the help of the computer random number generator. Remove new 
species from the old species. 

(9) step (i) from above 

(10) step (ii) from above 

(11) step (hi) from above 

(12) step (iV) from above 

(13) step (V) but now Nr concerns all species. 

(14) step (Vi), where N(t) concerns all species 

(15) step (Vii) 

(16) goto (6) unless the end of the simulation 



5 Results and discussion 

Limiting the amount of energy resources in the ecosystem introduces strong 
selection between the species competing for the resources. In our case the se- 
lection takes place through limiting the number of the new-born offsprings in 
the ecosystem. Namely, in the model, their number cannot exceed the value 
£po P {NE{t) /Q)Nr, where Nr is the number of parents, female and male, from 
all species. Thus, the individuals who are better adapted to the environment 
(larger value of Q, Eq.Q) dominate other individuals because they produce more 
offsprings. In our model, the better adaptation results from the phenotype x 
closer to the minimum x m i n of the function V{x,i). We could expect that in 
a time independent surrounding the mutations decreasing the fitness of the or- 
ganism lead to its elimination and in result to the elimination of the affected 
genes. However, in a changing environment the species need an ability to react 
to new situations. The quality of genes in a fluctuating surrounding is changing, 
i.e., we cannot tell any more that some genes are bad and some are good. There 
will be preferred genes which adapt better to the changes. 

In Fig.|2|, we show the effect of the changing environment on the old pop- 
ulation (cir = 5, L — 16, L' — 5), which was aging in a time independent 
surrounding through 40000 generations. The changing enviroment is modelled 
with the help of the coefficient b(t) in the function V(x,t) (EqJ|). We have 
chosen 

b(t) = -- -cos(wt), (10) 

where oj denotes the frequency of the environment changes. We should take 
into account that although we have a single species, in the example presented in 
FigM there are also present energy resources in the ecosystem which regenerate 
themselves according to EqJIl Thus, we have a kind of the Lotka-Volterra system 
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Figure 2: The effect of the transition from the time independent surrounding 
(&(*) = -1.17, t < 40000) to the one varying in time (b(t) = -\ - cos(O.OOlt), 
t > 40000) on the size of the single species (or = 5, a max = 11) and on the 
number of the energy units available for the species. The bottoom curve follows 
the size changes of the species whereas the upper curve shows the corresponding 
energy changes in the ecosystem. 



with the species as the predator and the energy as the prey. Once the fluctuating 
environment changes fitness Q of all individuals, the number of the new-born 
offsprings follows the environment fluctuations and in consequence there appear 
the energy oscillations. Notice, that in a changing environment the size of the 
population under consideration undergoes much stronger fluctuations than in 
the time-independent case. It is connected with the changes of the inherited 
genetic information. 

In Fig.|| we have presented the age profiles of the fraction of the inactive 
"chronological" genes in the population from Fig.||. The profiles concerning the 
generations at t = 20000 and t = 40000 refer to time independent surrounding. 
Although we have introduced in the model the back mutations for the inactive 
genes (the bit-strings corresponding to inactive genes may mutate into another 
bit-string with probability p — 0.1) we can observe that, while time is increas- 
ing, the shape of the profiles tends to the one which is time independent. It 
determines the average life span of individuals which is ranging from the age 
a = 1 to a = 7. On the other hand, the profiles corresponding to time dependent 
environment show both shrinking and enlarging of the life span. Effectively, the 
population becomes much younger than during its previous evolution in the time 
independent surrounding. More frequent fluctuations of the surrounding make 
this effect even stronger. It is evident from Fig.|] where the same old species, 
under consideration, experiences the surrounding changing with the higher fre- 
quency u = 0.01. We can observe that the effect of the frequently changing 
environment becomes more destructive for the older members of the species. 
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genotype site specific for age a 

Figure 3: The effect of the change of the time independent surrounding to 
the one varying in time, at t > 40000, on the distribution of the defected 
"chronological" genes for the situation presented in Fig.pl 
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Figure 4: The same as in in Fig.^ but u> = 0.01. 
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Figure 5: The phenotype evolution for the situation presented in Fig j^. 

All our simulation were done with a constant mutation rate (one mutation per 
genome) and we did not discuss the effect of the various values of the mutation 
rate on the above results. 

The evolution of the average phenotype of the species from the example in 
Fig.|| has been shown in Fig.||. The same data, but presented with the help of 
a histogram have been shown in Fig|| 

Notice, that there is a trend to minimize the environment fluctuations and 
the genes, which were well adapted in the time independent environment, do 
not have the ability to follow the environment fluctuations. They can be even 
eliminated from the genetic bank of the species if the fluctuations become more 
frequent. 

The environmental changes influence the size of all species in the ecosystem 
end this is the reason that some of them may be eliminated. The regeneration 
of the energy resources in the ecosystem takes some time and too rapid a growth 
of some species, e.g., mutant species, acts as a sudden cataclysm for the other 
species. The Lotka-Volterra systems have a self-regulatory character and there 
exist threshold values for the fraction of destroyed population above which the 
system returns to its previous state (]l9J,[p[). However, in the case of genetic 
populations, if the catastrophe is applied for many generations, one can observe 
that the age profile of the fraction of defective genes in the population may loose 
its stability and next the species becomes extincted Thus the appearance 

of the new species, even for a relatively short period, can eliminate the old 
species. 

In our computer simulations we associate a new value, of the reproduc- 
tion age with every new species appearing in the course of evolution. The value 
is determined with the help of a computer random number generator. Thus, 
we have also the possibility to investigate the effect of the reproduction age on 
the adaptation of the species to the fluctuating surrounding. We have observed 
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Figure 6: The histogram of the average usage of the phenotype x in the evolv- 
ing species (from Fig.^) when there is time independent surrounding (r.h.m 
histogram), and when there is fluctuationg environment. 

that the events, representing the appearance of the new species with the small 
value an (an ~ 1), are usually represented by short duration "bursts" of the 
population size and they vanish as rapidly as they have appeared. In the exam- 
ple from Fig.0 the are represented by the highest blobs painted in black, at the 
bottoom of the figure. It is easy to explain the phenomenon, as in this case the 
life span of individuals practically is shrinking to the activity of single gene and 
the individuals die after they reproduce themselves. The active genes cannot 
effectively adapt to the changing environment. How important the range of life 
span for species stability is, can be also concluded from the histograms in Fig.||. 

There are presented histograms of the average usage of the phenotype x of 
the species representing the predecessor (an = 5) and two descendant species 
(an = 2,7). It is evident that the average value x representing the species 
with an = 2 oscillates far away from the optimum values (variable minima 
of V(x,t)), because there is an insufficient number of genes adapted to the 
variable surrounding. The situation is a little bit better with the predecessor 
(old species). However, it is evident that it cannot be stable over a long time 
in a variable surrounding. The most stable species is the one for which an = 7, 
because in its population there are inherited genes which have adapted to the 
variable surrounding (the range of a; is almost symmetric with respect to x = 0). 

In this study, we did not discuss the mechanism of speciation in the new 
species. We have restricted our analysis to the "first order" speciation originat- 
ing from old species and we discussed in detail its effect on the inherited genetic 
information. We could expect a chain of events representing speciation with 
the species being better and better adapted to the environment. The younger 
the species the more probable speciation, and one should observe peaks in the 
number of the speciation events in a short period. 



12 




time (t) 



Figure 7: The effect of the speciation from the old species (earlier aging for 
40000 generations in time independent environment) in a fluctuating environ- 
ment (b(t) = — h — cos(O.Olt) on the evolution of the energy resources. In 
the bottom part of the figure there have been shown a few examples of the new 
species created in the course of the evolution. The energy cost of the new species 
creation or extinction is reflected by the terrace-like changes in the upper curve. 
The survived species have not been shown for clarity of the picture. 
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Figure 8: Histograms of the phenotype usage in the environment changing with 
the frequency uj = 0.02 in the case of three species: the predecessor population 
(an = 5), and two descendant populations (an — 7, an = 2). The frequency of 
the environment changes, u> = 0.02. 
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6 Conclusions 



We have discussed a model of species evolution in an ecosystem where the en- 
ergy resources regenerate themselves according to a logistic differential equation. 
The model belongs to the class of the Lotka-Volterra systems where the energy 
resources represent prey and the species represent predator. The number of 
the species present in the ecosystem under consideration results from evolution, 
which is understoood as the interplay of the two processes only, mutation and 
selection. We observed that after long periods of 'stasis' there happen bursts of 
the new species which are able to live for thousands of generations. We have 
observed that in a variable surrounding the species, for which the reproduction 
age an is too small, they are very unstable even if their size substantially exceeds 
the size of the populations specific for other species. In our model, they usually 
cause the elimination of other species. Simultaneously, the resulting increase in 
the energy resources makes the next speciations possible. 

There are two approaches to the description of species evolution: the one 
with continuous evolutionary changes and the theory of punctated equilibrium 
[p0[ according to which the evolutionary activity occurs in bursts. It is 
often the case that the mathematical models concerning the species evolution 
are restricted to pure species dynamics consideration or pure genetics evolution. 
We have shown that the inclusion of the genetic information into the population 



dynamics makes the use of Bak-Sneppen [|20|, 21 extremal dynamics (SOC) 
possible. 
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